function [U,V] = bpmf(U,V,max_epoch, dim, filter, train_vec, test_vec)

N = dim.N;
M = dim.M;
D = dim.D;

f_u = filter.f_u;
f_v = filter.f_v;

mu0 = zeros(D,1);
b0 = 1;
W0 = eye(D);
nu0 = D;
alpha = 2;

fprintf('Gibbs sampler starting ...\n');

for epoch =  1: max_epoch
    tr_rmse = calc_rmse(U,V,train_vec);
    te_rmse = calc_rmse(U,V,test_vec);
    fprintf('epoch %2d %2.4f %2.4f\n',epoch,tr_rmse,te_rmse);
    
    % sample user hyperparams
    U_bar = mean(U,2);
    mu0_ = (b0*mu0 + N*U_bar)/(b0+N);
    b0_ = b0 + N;
    nu0_ = nu0 + N;
    W0_ = inv(inv(W0) + U*U' - N*(U_bar*U_bar')  + ...
        (b0*N/(b0+N))*(mu0-U_bar)*(mu0-U_bar)');
    W0_ = (W0_+W0_')/2;

    L_U = wishrnd(W0_,nu0_);
    iW0_ = inv(b0_*W0_);
    mu_U = mvn_rnd(mu0_,iW0_);
    
    % sample movie hyperparams
    V_bar = mean(V,2);
    mu0_ = (b0*mu0+M*V_bar)/(b0+M);
    b0_ = b0 + M;
    nu0_ = nu0 + M;
    W0_ = inv(inv(W0) + V*V' - M*(V_bar*V_bar') + ...
        ((b0+M)/(b0+M))*(mu0-V_bar)*(mu0-V_bar)');
    W0_ =(W0_+W0_')/2;
    
    L_V = wishrnd(W0_,nu0_);
    iW0_ = inv(b0_*W0_);
    mu_V = mvn_rnd(mu0_,iW0_);
    
    for ii = 1:N
        jj = f_u{ii}.jj;
        val = f_u{ii}.val;
        if ~isempty(jj)
            VV = V(:,jj)*V(:,jj)';
            VR = V(:,jj)*val;
            L_i = L_U + alpha*VV;
            Li_i = inv(L_i);
            mu_i = Li_i*(L_U*mu_U + alpha*VR);
            U(:,ii) = mvn_rnd(mu_i,Li_i);
        end
    end
    
    for jj = 1:M
        ii = f_v{jj}.ii;
        val = f_v{jj}.val;
        if ~isempty(ii)
            UU = U(:,ii)*U(:,ii)';
            UR = U(:,ii)*val;
            L_i = L_V + alpha*UU;
            Li_i = inv(L_i);
            mu_i = Li_i*(L_V*mu_V + alpha*UR);
            V(:,jj) = mvn_rnd(mu_i,Li_i);
        end
    end
end